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1  Introduction  and  Overview. 


This  work  investigates  the  attachment  landing  problem  for  flapping  wing  micro-aerial  vehicles 
(MAV’s).  Landing  is  a  critical  operational  component,  necessary  for  conducting  surveillance,  and 
conserving/harvesting  energy.  Autonomous  landing  requirements  for  the  MAV  are  quite  different 
than  for  conventional  aircraft  as  the  MAV  may  be  expected  to  identify  its  own  set-down  point 
while  operating  in  highly  constrained  and  unfamiliar  surroundings.  Recent  advancements  in  surface 
adhesion  technology  offer  to  greatly  enhance  landing  capabilities  by  allowing  the  MAV  to  attach  to 
any  nearby  vertical  surface,  much  like  flying  insects,  and  then  detach  with  little  effort.  Researchers 
have  recently  attempted  to  design  (Jagota  and  Bennisson,  2002;  Glassmaker  et  al.,  2004;  Hui  et 
al.,  2004;  Bhushan  et  al.,  2006)  and  fabricate  (Geim  et  al.,  2003;  Sitti,  2003;  Sitti  and  Fearing, 
2003;  Northen  and  Turner,  2005;  Yurdumakan,  et  al.,  2005)  these  types  of  adhesives,  also  called 
bio-inspired  tapes. 

Biological  studies  have  revealed  that  flying  insects  execute  landing  maneuvers  ranging  in  com¬ 
plexity.  For  example,  the  honeybee  utilizes  a  very  simple  control  strategy  when  landing  upright  on 
a  flat  horizontal  surface  (Srinivasan  et  al.,  2001).  Alternatively,  the  housefly  will  land  upside-down 
on  a  ceiling  by  extending  its  forward  legs  over  its  head,  making  contact,  and  then  using  momentum 
to  rotate  the  remainder  of  its  body  180-degrees  to  the  ceiling  (Borst,  1990).  The  particular  landing 
strategy  adopted  by  insects  is  dependent  on  various  factors  including  inherent  flight  capabilities  and 
“landing  gear”;  the  honeybee  is  much  more  adept  at  hovering  flight  while  the  housefly  has  superior 
surface  adhesion.  The  agility  demonstrated  by  flying  insects  is  difficult  to  mimic.  The  equations 
of  motion  for  the  flapping  wing  vehicle  are  naturally  unstable  (a  key  factor  of  agility)  and  contain 
substantial  nonlinearities.  The  control  problem  has  been  recently  investigated  by  relatively  small 
number  of  researchers  (Deng  et  al.,  2006;  Schenato,  2003;  Schenato  et  al.,  2003).  The  standard 
approach  to  control  is  based  on  time  averaging  theory  (Sanders  and  Verhulst,  1985),  which  stems 
from  the  fact  that  the  motion  of  the  wings  is  substantially  faster  than  the  motion  of  the  body  to 
which  they  are  attached. 

The  problem  addressed  by  this  work  that  of  vertical  attachment  landing  where,  generally  stated, 
the  objective  is  to  control  the  vehicle  such  that  its  attachment  surface  can  make  contact  with  a 
vertical  surface.  It  is  assumed  that  the  contact  surface  is  situated  on  the  underbelly  of  the  vehicle, 
so  that  the  landing  requires  a  relatively  large  pitch  angle.  The  two  primary  questions  addressed 
here  are:  (1)  what  are  the  actuation  requirements  to  produce  the  required  landing  maneuver?;  and 
(2)  how  can  this  maneuver  be  controlled? 

The  outline  of  this  report  is  as  follows.  The  equation  of  motion  for  a  flapping  wing  vehicle  are 
derived  in  Sec.  2,  with  the  aim  of  generality.  In  Sec.  3,  the  longitudinal  equations  are  considered 
with  the  objective  of  deriving  general  expressions  for  the  actuation  requirements  of  the  vehicle.  In 
this  section,  several  critical  assumptions  as  to  the  necessary  wing  degrees  of  freedom.  Methods  of 
nonlinear  control  are  applied  to  the  longitudinal  dynamics  in  Sec.  4.  Section  5  contains  a  numerical 
example  of  the  theoretical  development  of  the  preceding  sections,  based  on  the  planform  properties 
of  the  bumblebee.  Conclusions  of  this  work  are  stated  in  Sec.  6. 

2  Flight  mechanics. 

2.1  Kinematics. 

The  body  is  divided  into  three  components:  the  main  body,  and  the  two  wings.  The  main  body 
kinematics  are  referenced  to  an  inertial  reference  frame,  denoted  ^  required  by  the  kinetic 
equations.  The  motion  of  the  wings  are  referenced  to  the  main  body  reference  frame,  ^i.  The 
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notation  followed  throughout  this  section  is  explained  in  the  Appendix  (Sec.  6.1). 

2.1.1  Main-body  kinematics. 

Let  the  position  vector  R  locate  a  material  point  of  the  flight  vehicle  relative  to  an  inertial  reference 
frame;  [R]o  is  the  ,^o- represent  at  ion  of  this  position  vector.  Let  Ri  locate  the  origin  of  a  reference 
frame  that  is  attached  to  some  rigid  portion  of  the  main  body  and  confined  to  rotate  and 
translate  with  it.  A  material  point  is  located  relative  to  by  the  vector  R'  so  that 


R  =  Ri  R  . 


(2.1) 


The  ^0  frame  representation  of  this  relation  is 


[-R]o  -  [i^ilo  +  =  [-Ri]o  +  Coi[R']i, 


(2.2) 


and  time  derivative  is 

A([i2]o)  =  Dt{[Ri]o)-Coi[RTi[^i/o]i+CoiDt{[R']i).  (2.3) 

The  second  time  derivative  is 

At([^]o)  =  ^tt{[Ri]o)  —  Coi[u)ifQ]l[R']l[iVi/Q]i  —  CQi[R']lDt{[iVi/o]i) 

-f  2Coi[u;i/o]I A([i^']i)  -f  CoiDaimi).  (2.4) 

Define  vi  =  CioDt{[Ri]o)  so  that 

At([^i]o)  =  C'oiA(^i)  -b  Coi[ui/q]\vi.  (2.5) 

Also,  let  vi  =  Dt{vi),  uji  =  [u;i/o]i,  a>i  =  [^i/o]v  =  A(i^i),  =  A([i^']i),  and 

R'  =  [R']l,  so  that 

CioDt{[R]o)  =  vi-R'u;i  +  R'  (2.6) 

CioDii{[R]o)  =  vi  QiVi  —  CjiR^lji  —R'lji  —2R^iVi  R\  (2*7) 

These  are  the  two  expressions  required  to  formulate  the  dynamic  equations.  They  are  expressed 
relative  to  a  fixed  body  frame.  However,  since  the  aerodynamic  forces  are  expressed  relative  to  the 
direction  of  air  flow,  it  is  often  convenient  to  express  the  dynamics  relative  to  a  frame  that  aligns 
with  the  air  flow. 

The  body  velocity  is  Vi  :=  dRi/dt^  hence  Z)i([ili]o)  =  CoiVi,  For  the  circumstance  that 
atmosphere  is  non-stationary  (i.e.,  the  presence  of  wind  gusts  and  cross-flow),  the  wind  velocity, 
denoted  is  necessary  for  an  accurate  prediction  of  the  aerodynamic  forces.  To  accommodate 
this  situation  the  flight  equations  can  be  constructed  in  terms  of  the  relative  velocity 


Vl/u;  =  Vi  -  (2.8) 

which  is  the  negative  of  the  total  wind  velocity.  The  component  Vyj  is  not  a  known  quantity  and 
thus  flight  control  designs  are  based  on  the  assumption  that  it  is  negligible.  The  ^o-representation 
of  the  previous  equation  is 

[’^l/«;]o  =  [Vi]o  -  [V^]o  =  A([i2l]o)  -  [V^]o  =  Com  -  [V„]o. 
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(2.9) 


If  the  disturbance  term  Viy  is  zero,  then  it  is  clear  that  vi  =  C'lofV'i/^Jo  = 

The  relative  wind-axis^  denoted  is  defined  such  that 

C-n[Vyw]i  =  [VT,  0  Of, 

where  :=  and  C\i  is  the  basis  transformation  operator  ► 

The  relative  motion  of  and  results  in  the  angular  velocity  vector  which  describe 


the  relative  angular  motion  of  the  reference  frames. 

Taking  the  time  derivative  of  the  relative  velocity  vector,  results  in 

A([^l/ti;]o)  =  +  C'oi  l)-  (2.10) 

Since  A([^i/i/;]o)  =  At([A]o)  -  A([^u;]o),  it  follows  that 

v\  -  CioDt{[Vyj]o)  =  Dt{[Vi/^]i)  (2.11) 

It  is  noted  that 

A(Ai['^i/t.]i)  =  +  AiA([Vi/^]i).  (2.12) 

For  no  disturbance,  v  =  and  vi  =  A([^i/it;]i);  hence 

CiiVi  =  A(Al[^l/'Lt;]l)  +  [^l/l]lAl[^l/'u;]l-  (2.13) 

Since  Cii[Vi/^]i  —  [Vr^  0  0]^,  it  follows  that  A(C'ii[^i/tij]i)  =  [^t  0  0]^i  and  thus 

C-nVi  =  [VT,  0  Of +  0  Of.  (2.14) 

The  wind-axis  transformation  is  then 


Ai(^i+A^i)  =  [^1  0  of  ^[u;i/i]i[Vti  0  of  +  CiiclJivi 

=  [Vt,  0  of  +  [wi/ilf 0  Of +  C'iiwiC'ii[l^Ti  0  0]^  (2.15) 

=  t’l  +  ([wi/i]j  +  CiiwCii) 


where  the  symbols  uj  and  Uf  are  evident. 

The  basis  transformation  is  parameterized  by  the  rotation  by  the  3-2-1  set  of  Euler 

angles  0i  :=  [(/>i  di  'ipif .  The  transformation  matrix  is 


Ao  = 


c{6i)c{tpi)  c{ei)s{tpi)  -s{0\) 

-c(</)i)s(V’i)  +  s{4>i)s{ei)c{tpi)  c{(f)i)c{ipi)  +  s{(f)i)s(9i)s(rpi)  s(0i)c(6'i) 

s{4)i)s{ipi)  +  c{4)i)s{ei)c{ipi)  -s{4)i)c{tPi)  +  c{4>i)s{ei)s{ij;i)  c{(f)i)c{di) 


(2.16) 


The  relationship  between  cvi  and  0i  is  then  given  by 


UJl 


1  0  —  sin^i 

0  cos  (f>i  sin  (pi  cos  6i 
0  —  sinc^i  cosc^icos^i 


Bi. 


(2.17) 


The  relative  motion  between  and  ^i  is  quantified  by  the  3-2-1  set  of  Euler  angles  71  =  0,  ai 
(angle  of  attack),  and  A  (sideslip  angle),  in  which  case  the  transformation  matrix  is 


cos  A 

sin/Ji 

0  ■ 

cosai 

0 

sinai 

Cii  — 

—  sin  Pi 

cos  Pi 

0 

0 

1 

0 

0 

0 

1 

—  sinai 

0 

cosai 
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Let  Fi  =  [71  ai  pi].  The  components  of  the  representation  of  the  relative  angular  velocity 
vector  cjj/i  are  related  to  the  relative  angles  by 


0  —  sin  Pi  0 

0  —  cos  Pi  0 

0  0  1 


ri. 


(2.19) 


2.1.2  Wing  kinematics. 


The  wings,  considered  rigid,  are  attached  to  the  main  body.  To  quantify  the  motion  of  the  right 
wing,  two  fixed  reference  frames  are  defined,  ^2  and  ^3  (cf.  Fig.  1.).  The  origin  of  ^2  is  located 
at  a  wing-body  connection  point.  The  origin  of  is  located  at  the  wings  aerodynamic  center 
of  pressure.  These  two  frames  may  be  such  that  C23  ~  0,  but  C23  ^  /;  that  is,  there  is  some 
initial  offset  between  the  reference  frames.  The  X2  and  0:3  axes  are  directed  along  the  chord  of  the 
wing,  the  2/2  and  2/3  axes  are  directed  along  the  span  toward  the  right,  and  the  22  and  23  axes  are 
tangent  to  the  chord  and  span  directions  (following  the  right-hand  rule).  Let  V2  and  denote 
the  velocity  at  the  origin  of  ^2  and  <^3,  respectively.  For  the  left,  wing  the  frames  ^4  and  ,^5  are 
defined  in  an  identical  manner.  Generally,  Vi  denotes  the  velocity  of  the  origin  of  frame  ^i. 

The  orthogonal  transformation  matrices  C21  and  C41  quantify  the  orientation  of  ^2  and  ^4, 
respectively,  relative  to  ^1.  Again  utilizing  the  3-2-1  set  of  Euler  angles,  the  transformation 
matrices  are  given  by 


Cii  = 


c{dl)c{xpi) 

-c{(i>'i)s{i>'i)  +  s{(})'ps{e'pc{i)'p 
s{<}>'ps{xp')  +  ci4>'psie')c{i^') 


for  i  =  2,4,  where 


c{di)s{rpi) 

+  s{4>'-)s{e'ps{tjj[) 


s{4>dc{0') 

s{4>'M'>Pi)  +  c{4>'i)s{0'i)s{tp'i)  c{<j)'pc(e'p 


4>'i 

0'i 

^'i 


(f>i  +  4>i 

0^  +  0i 


(2.20) 


Here  4>^,  Of,  tpf  are  constants  quantifying  the  initial  offset.  The  angles  (^>2  and  (f>3  are  called  the 
flapping  angles.  The  angles  The  angles  O2  and  0^  are  called  the  pronation  angles.  Finally,  'tp2  and 
tpi  are  called  the  stroke  angles. 

Let  (f>3,  63,  xp3  denote  the  constant  offset  parameters  of  C32.  The  resulting  transformation  is 


<^32 


c(03)c(V>3)  C{e3)s{lp3)  -siOs) 

-c{(f>3)s{rp3)  +  s{(f>3)s{03)ci'tp3)  c{(f)3)c{lp3)  +  s{(f>3)s{03)s{tp3)  s{(f>3)c{03) 

s{(f>3)s{'(p3)  +  c{(j)3)s{03)c{rp3)  -s{<f)3)c{ll>3)  +  c{(f>3)s{03)s{lp3)  c(<^3)c(03) 


Similarly,  for  <754: 


<^54  = 


c(^5)c(V’5) 

-C{4>5)s{tj)5)  +  s(<^5)s(05)c(l/>5) 
s{(f>5)s{lp5)  +  c(<f>5)s(05)c(l/>5) 


c(6’5)s(V’5)  -5(6*5) 

C{(f>5)c{tjj5)  +  s(<^5)s(05)s(t/’5)  s{(j)s)c{9s) 

-s{(j)3)c{i)s)  +  c{(i>s)s{es)s{'(l)s)  c(05)c(05) 


(2.21) 


(2.22) 


The  relative  rotation  rate  of  the  reference  frames  is  quantified  by  the  vectors  u;2/i  and  ^4/1.  Let 
^i/\  =  and  ©i/i  =  [<f)\  e[  for  i  =  2,4  so  that 


a>t/i  = 


1  0  —  sin  0[ 

0  cos  4>\  sin  (/)  ■  cos  6\ 
0  —  sin  (/)(  cos  cos  0^ 


9,; 


t/i 


G,:/,0 


(2.23) 
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Define  two  additional  reference  frames  ^3  and  <^5,  which  indicate  the  relative  orientation  of  the 
wind  respect  to  the  motion  of  the  wings.  The  origins  of  these  frames  are  coincident  with  <^3  and 
^5. 

The  wing  velocity  relative  to  the  wind  is  given  by 


=  Vi 


3,5 


where  V if is  the  velocity  of  relative  to  the  airflow.  Define  (z  =  3, 5)  such  that 


Qi[^i/Ji  =  0  of ,  i  =  3,5 


(2.24) 


(2.25) 


where  C33  and  are,  respectively,  the  basis  transformation  operators  ^3  — +  ^3  and  ^5  — >  <#5. 
These  two  transformations  are  parameterized  by  3-2-1  Euler  angles  7^  =  0,  and  j3i;  the  latter  two 
angles  are  the  wing  angle  of  attack  and  sideslip  angle.  The  orthonormal  transformations  matrices 
are 


Cii  = 


—  sin  a[  cos  —  sin  a'  sin  pi  cos 

—  sin  Pi  cos  Pi  0 

—  cos  a[  cos  P  cos  sin  P  —  sin  j 


(2.26) 


where  =  ‘Kjl  —  a.  The  relative  rotations  of  ^i  and  ^i  are  quantified  by  the  angular  velocity 
vectors  u?3/3  and  0^5/5,  with  the  relation 


0  —  sin  Pi  0 

0  —  cos  pi  0 

0  0  1 


(2.27) 


where  F*  =  [7^  ai  i  =  3, 5  . 

2.2  Kinetics. 

The  equations  of  motion  of  an  arbitrary  nonpolar  deformable  body  can  be  expressed  in  the  general 
form 


F*  -F-F'  =  0, 


(2.28) 


where  F*  the  generalized  inertial  force^  F  the  generalized  applied  force^  and  F'  the  generalized 
constraint  force.  Explicitly,  these  terms  are 


F*  = 

(2.29) 

F  = 

(2.30) 

F'  = 

(2.31) 

where  is  an  inertial  frame,  y  is  an  array  of  coordinates  chosen  to  specify  the  motion  of  the  body, 
V  is  the  spatial  volume  of  the  body,  and  b  and  b'  are,  respectively,  the  applied  and  constraint  force 
per  unit  mass.  Constraint  forces  do  not  typically  show  up  in  flight  analysis  since  the  aircraft  is  a 
free  body.  However,  they  can  be  used  to  determine  the  internal  actuation  forces  due  to  specified 
motions  (i.e.,  morphing)  of  the  aircraft. 
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When  the  main  body  is  rigid,  the  kinematic  parameters  v\  and  oJi  adequately  quantify  its  motion. 
These  kinematic  parameters  are  regarded  as  quasi- cordinates  since  their  integration  does  not  yield 
any  physically  meaningful  quantity.  They  are  related  to  underlying  kinematic  quantities  by 

vi  =  CioDt{[Ri]o) 

CJi  =  Gi©!. 


If  the  body  is  not  rigid,  additional  coordinates  are  required  to  specify  its  motion.  These  additional 
coordinates  can  be  treated  in  one  of  two  ways:  as  additional  states  of  the  dynamical  equations,  or 
time-dependent  extraneous  inputs.  For  the  analytical  evaluation  of  complex  structures,  the  latter 
approach  is  typically  more  useful.  Mathematically,  this  amounts  to  the  specification  of  program 
constraints,  or  servoconstraints,  on  certain  kinematic  parameters. 

Likewise,  the  translational  and  rotational  velocity  of  the  (assumed  rigid)  wings  are  sufficient  to 
describe  their  motion.  Additionally,  there  must  be  constraint  equations  that  couple  the  wing  to  the 
main  body.  However,  a  significant  simplification  ensues  when  the  wings  are  considered  massless. 
In  this  case,  the  generalized  inertial  force  is  dependent  solely  on  main  body  motion,  and  also  the 
generalized  constraint  force  goes  to  zero.  Defining  the  coordinate  array 

y=[vi  wi  ,  (2.32) 


the  generalized  inertial  force  can  be  written 


Fj*  =  m{v  -\-  uv)  —  Suj  —  LuSuj  —  25u;  +  S 
^2  =  S  {v  Qv)  +  Juj  +  CjJuj  -h  Jo;  +  intj)pp^  dm 


where  p  =  [R^]i,  F*  =  [F*  F^]'^,  and 


m 

S 

J 


L 

L 


dm 

p  dm  =  fc 


pp  dm 


(2.33) 

(2.34) 


are,  respectively,  the  mass  and  the  first  and  second  moments  of  inertia. 

The  matrix  form  {i.e.,  reference  frame  dependent)  of  the  translational  and  rotational  equations 
of  a  morphing  aircraft  can  be  expressed  in  the  form  (Seigler  et  al.,  2007) 


mv  +  muv  —  mrc^  —  ujmfco)  —  2mfcaj  rc  —  F  (2.35) 

mfcv  +  mfcdjv  +  Juj  +  ljJuj  -h  Juj  +  int^PP  dm  =  A4  (2.36) 

where  the  components  of  the  generalized  applied  force  are  F  =  [F  The  equations  reduce  to 

the  standard  rigid  body  flight  equations  when  p  =  0  for  every  material  point  and  the  center  of  mass 
is  chosen  as  the  reference  origin,  requiring  that  rc  =  0.  For  the  non-rigid  case,  it  is  also  possible 
to  let  the  body- fixed  axis  move  with  the  center  of  mass,  in  which  case  again  rc  =  0.  However,  v 
would  then  be  the  velocity  of  the  center  of  mass  rather  than  a  “fixed”  point  of  the  vehicle.  For 
sufficiently  large  aircraft  speeds,  this  diflFerence  will  have  only  a  small  effect.  Note  also  that  the 
moments  of  inertia  may  become  more  complicated  when  the  reference  point  is  in  relative  motion. 
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For  the  specified  generalized  velocities,  the  generalized  applied  forces  are 

=  /  [b]i  dm=[f“]i  +  [f9]i 

Jt) 

M  =  [  p[b]i  dm  =  [m“]i  +  [m®]i, 

Jv 


(2.37) 

(2.38) 


which  are  divided  into  aerodynamic  and  gravitational  components.  The  gravitational  force  is  given 

by 

[f®]i  =  Cio[0  0  mgf. 


The  moment  due  to  gravity  is 

[m®]i  =  Ciorc[0  0  mgf . 

A  direct  conversion  of  the  inertial  force  to  the  wind  axis  results  in 

C\iF*  =  TTi  I'Uj  +  —  C\\Suj  —  C\iIjjSuj 

—  “h  C\\  ’S 

and 

C\\F2  =  S  +  C\iu)Ciyj)  vi}  +  CiiJuj  +  Cii  (^ujJ  +  ( 

+  Chintvpp  dV 


(2.39) 

(2.40) 

(2.41) 

(2.42) 


Thus,  ignoring  the  presence  of  wind  gusts  and  cross- flow  etc.,  the  equations  of  motion  in  the  wind- 
axis  become 

mvi  -1-  m  (a>i/i  +  ChCbCii)  uj  —  -  mCiiujrc^  —  2mCjifca;  +  mC\ifc  = 

mrcvi  +  mfc  +  CnJCo  +  Cn  {uJ  lj  Cii  J  pp  dV  =  C\iM. 

2.3  Aerodynamic  Forces 

The  aerodynamic  force  applied  to  a  wing  is  quantified  by 


ff=  /  tdSu  i  =  3,5 

is^ 


(2.43) 


where  is  resultant  aerodynamic  force  applied  at  the  origin  of  located  at  the  center  of  pressure, 
t  is  the  traction  force  per  unit  area,  and  Si  denotes  the  external  area  of  the  wing.  The  moment  of 
the  aerodynamic  force,  about  the  center  of  mass  of  the  vehicle,  is 


m"  =  J  Ri^i  X  t  dS, 


(2.44) 


where  Ri/c  locates  the  origin  of  located  at  the  center  of  pressure  of  the  wing,  with  respect  to 
<^1. 

The  aerodynamic  forces  are  considered  to  be  functionally  dependent  on  the  instantaneous  position 
of  the  wing  and  its  velocity,  both  relative  to  the  airflow.  The  aerodynamic  force  vector  is  divided 
into  three  components 


r  =  ff  -h  +  f5^ 


(2.45) 
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which  denote,  respectively,  the  force  on  the  main  body,  the  right  wing,  and  the  left  wing  (the  indices 
correspond  to  the  reference  frame  where  the  resultant  force  is  applied).  Further  denote 


[ffli  =  [AOLif  (2.46) 

[f3“]2  =  [DsOLz]^  (2.47) 

[f5“]g  =  [DiOL^f  (2.48) 

where  L  and  D  are  the  lift  and  drag  forces.  If  it  is  reasonable  to  assume  that  the  aerodynamic 
forces  on  the  main  body  are  small  compared  to  the  wings,  which  is  the  case  for  “slow  flight”,  it  can 
be  assumed  that  is  negligible.  The  functional  form  of  the  lift  and  drag  is  expressed  as 

Li  =  Li{au  OLu  Pu  Pu  Pi,Qu  V'^^)  (2.49) 

Di  =  Di{oii^6ii,  pi^  Pi^  Pi^  Qi^V^*).  (2.50) 


for  i  =  3,5,  where  Pi  and  Qi  is  a  component  of  [u^i/o]i  =  [Pi  Qi  Ri]'^ • 

Based  on  the  published  literature  (Sane,  2003;  Sane  and  Dickinson,  2001,  2002;  Dickinson  et  al., 
1993)  regarding  the  aerodynamics  of  flapping  wings,  the  instantaneous  aerodynamic  forces  applied 
to  the  wings  are  composed  of  three  distinct  phenomena:  delayed  stalls  rotational  lift,  and  wake 
capture.  Delayed  stall  contributes  a  significant  majority  of  the  aerodynamic  force.  Rotational  lift 
and  wake  capture  are  due  to  the  pronation  of  the  wing,  which  typically  occurs  at  the  peaks  of 
the  up  and  down-stroke.  Subsequent  development  is  based  on  the  assumption  that  wake  capture 
is  the  predominant  force,  while  rotational  lift  and  wake  capture  are  negligible.  The  wake  capture 
component  is  represented  by  a  force  located  at  the  center  of  pressure  of  the  wing  that  is  dependent 
on  the  instantaneous  angle  of  attack,  q^,  sideslip  angle,  pi,  and  the  relative  wind  velocity,  Vj-.: 

it  =  if{VT,,ai,pi).  (2.51) 

The  corresponding  moment  is  then 

mf  =  Ry,xif{VT,,ai,Pi). 

Moreover,  the  force  and  moment  are  defined  in  the  local  wind  axis;  i.e., 

Kh  =  [ii,/1  xf“(W„Qi,A)h 

=  {Ri/i]imVT,,ai,mi-  (2.52) 

In  the  frame,  the  moment  is 

K]i  =  CHQ[i^,/l]^ff(FT„ai,ft)k.  (2.53) 

The  position  vector  for  the  right  wing  is  given  by 

[^3/1)3  =  [^2/1)3  +  [-^3/2)3 

=  C'33C32C2i[H2/i]i  +  C'33C32[H3/2]2 

where  [i22/i]l  [■^3/2)2  are  constants.  Similarly,  the  position  of  the  left  wing  is  given  by 

[■^s/lls  =  C'55C'4i[^4/i]i  +  ^^55  054(^5/4)4  (2.54) 

where  [i^4/i]i  and  [^5/4)4  are  constants. 
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The  aerodynamic  force  is  dependent  on  the  velocities  of  the  wings,  V3  and  V5,  which  are  given 


by 

[V’aja  =  <^31  +^1  ([i?2/l]l  +  C'i2[-R3/2]2)]  +  C'32[‘*’2/i]2[-^3/2]2  (2.55) 

[V’sjs  =  C'si  [^1  +t^r  ([-R4/i]i  +  +  C'st [‘^4/1)4 (2.56) 

When  the  main  body  is  stationary,  or  else  if  the  first  term  of  the  summation  is  much  less  than  the 

second  (i.e.,  the  body  velocity  is  slow  relative  to  the  wing  velocity),  this  reduces  to 

[V'aja  =  C32[tU2/i]2[i23/2]2  (2.57) 

=  C'54[‘<^4/l]4[-^5/4]4-  (2.58) 

The  velocities  on  which  the  aerodynamic  force  depends  is  now  determined  by 

Vpi 

0  =CiAyi\u  1  =  3,5.  (2.59) 

0 


Assuming  that  [Viji  is  known,  this  constitutes  three  equations  with  three  unknowns:  Qj,  /3i,  V^.. 


3  Actuation  requirements. 

In  this  section  are  methods  of  determining  the  actuation  requirements  of  attachment  landing, 
without  consideration  of  feedback  control.  Also  introduced  in  this  section  are  several  key  simplifying 
assumptions  that  are  necessary  to  achieve  analytically  tractable  results. 

3.1  Aerodynamic  forces  in  “slow^”  flight. 

When  the  flight  speed  of  the  vehicle  is  “slow”,  it  is  reasonable  to  assume  that  v\  =  uj\  =  0;  and 
thus  Eq.  (2.59)  reduces  to 


■  VV3  ■ 

sin  03  cos  /?3  sin  Oz  +  cos  03  cos  Oz 

0 

—  ^33^32  [^2/1)2  [-^3/212  — 

sin  Pz  sin  ^3 

0 

_  cos  03  cos  (3z  sin  ^3  —  sin  03  cos  Oz 

It  follows  that 

03  =  ^3  +  7r/2 
133  =  0 

=  2/3/202- 

Similarly,  for  the  left  wing 

05  =  ^5  +  7r/2 

/?5  =  0 

=  2/5/404- 

According  to  Dickenson  et.  al  (1999),  the  magnitude  of  the  aerodynamic  force  can  be  accurately 
quantified  by 

Li^Vxi^cki)  =  —pSiV'p^Ci  sin  i  =  3,5,  (^-1) 
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where  ai  G  [0,7r/2],  and  Ci  is  a  positive  constant  (a  value  of  3.5  was  determined  empirically  by 
Dickenson  et.  al,  1999).  For  “slow”  flight,  this  results  in 

LziVr^.as)  =  -pSs{yy2^2)^Css\n{9s  +  7r/2)  (3.2) 

=  2P‘55(y5/404)^C5sin(05  +  7r/2).  (3.3) 

Note  that  the  lift  is  in  the  direction  opposite  of  the  wing  motion.  Expressed  in  the  frame,  the 
aerodynamic  forces  are 

[f"]l  =  C'i3[f3  ]3  +  Cislfsls 

=  Ci2C23C35[f3"]3  +  Ci4C45C55[f5"]5  (3.4) 

The  forces  are  dependent  solely  on  the  motion  of  the  wings  relative  to  the  body. 

The  aerodynamic  moment  due  to  the  force  is  found  by 

Kll  =  [«3/l]I[f3ll,  (3.5) 

where 

[■^3/l]l  =  [■^2/l]l  +  C'i2[.R3/2]2.  (3.6) 


3.2  Method  of  averaging. 

If  the  forces  and  moments  produced  by  the  wings  occur  periodically,  and  at  a  time  scale  much  less 
than  the  motion  of  the  vehicle  body,  the  method  of  averaging  can  be  used  to  produce  a  substantially 
simplified  dynamic  model.  To  show  this,  consider  a  periodic  dynamical  system  of  the  form 

x  =  f{x,t),  (3.7) 

where 


f{x,t)  =  f{x,t  +  T) 


for  some  constant  T  >  0.  The  corresponding  average  systems  is 

cT 


1 

^-avgit)  ~  ^  J 


fdt. 


With  a  change  of  timescale,  i.e.,  t  =  rT,  Eq.  (X)  can  be  written 


|  =  r/(x,.r), 


and  thus 


dx, 


avg 


dr 


Tf{^avg)- 


(3.8) 


(3.9) 

(3.10) 


It  is  desired  to  approximate  the  dynamical  system  of  Eq.  (3.13)  with  that  of  the  corresponding 
averaged  system,  Eq.  (3.14).  The  theory  of  averaging  (see  e.g.,  Sec.  10.4  of  Khalil,  2002)  states 
that  if  T  is  sufficiently  small  then  x{t)  —  Xavg{t)  =  0(T),  where  O  denotes  the  order  of  the  error. 
Hence,  the  difference  between  the  output  of  a  dynamical  system  and  its  averaged  output  is  small 
when  the  periodic  frequency  is  large.  This  result  can  be  extended  to  the  case  where  the  system 
contains  a  state- dependent  control  input;  i.e.,  x  =  /(x,  u(x),t)  (Schenato  et.  al.,  2003). 

Since  the  wing  motion  of  a  flapping  wing  vehicle  is  generally  much  faster  than  the  gross  motion  of 
the  vehicle,  the  method  of  averaging  is  assumed  applicable.  All  subsequent  theoretical  development 
is  based  on  this  assumption.  Its  validity  will  be  examined  in  the  numerical  example  of  Sec.  5. 
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3.3  Average  forces  and  moments. 

The  wing  motion  is  separated  into  the  down-stroke  and  the  up-stroke.  The  down-stroke  is  assumed 
to  take  place  over  a  period  Ti  seconds,  and  the  up-stroke  over  a  period  of  T2  seconds;  the  total  period 
is  T  =  Ti  +  T2.  During  downstroke  and  upstroke,  the  02  traverses  the  same  angular  range,  at 
an  assumed  constant  angular  velocity.  The  downstroke  and  upstroke  angular  velocity  magnitudes 
are  denoted  02  ^2  (both  are  positive  v'alued).  The  periods,  amplitude,  and  angular  velocities 

are  related  by 


Let  denote  the  flapping  frequency,  which  is  related  by 


^(t>U  = 


(3.11) 


(3.12) 


On  the  down  stroke,  let  ^3  =  — 7r/4  so  that  =  7r/4;  it  follows  that 

^12^23^3  =  C'12 


-10  0 
0  1  0 
0  0-1 


The  force  on  the  down-stroke  is  thus 

1 


[fsli  =  —2p^3iy3/2<l>2)  C sin(03  +  7r/2) 

On  the  upstroke,  let  83  =  7r/4  so  that  a'^  =  it  follows  that 


s4>2S'1{)2  +  C(j)2s82Ctl)2 
-s4>2Clp2  +  C(j>2s92S'tp2 
c4>2Cd2 


(3.13) 


C12C23C33  =  C12 


10  0" 
0  1  0 
0  0  1 


The  force  on  the  up-stroke  is  thus 


[fall  =  ^pS3{y3/2^2fC  sin{e3  +  Tr/2) 


S<j)2S1p2  +  c4>2S92Clp2 

-S(l>2Clp2  +  C(j>2s92S'lp2 
c4>2c92 


(3.14) 


Focusing  on  the  down-stroke,  the  angle  02  is  taken  to  be  constant,  and  62  =  0.  Thus,  02  is  the 
only  non-constant. 

The  average  force  on  the  down-stroke  is 


S02  ^^2  dt  +  s92Ctp2  C02  dt 

-C02  ^^2  dt  +  S02S02  lo^  ^^2  dt 

^  062  C02  dt 

where 

\/2 

Q  =  —pS3{y3/2^)^C. 


(3.15) 
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Similarly,  the  average  up-stroke  force  is 


Q 

2T2 


Slp2  502  dt  -f  S62C1P2  Jt  ^^2  dt 

— C02  Jti  ^^2  dt  +  5^2502  Iti  ^^2  dt 
002  C02  dt 


The  average  force  over  a  cycle  is  found  by 

/ 


'2T 


VL 


50^  502  dt  +  s92Cfll^2  Ip  ^^2  dt 

-<^2  lo^  s4>2  dt  +  S^2.‘?V’2  lo^  c4>2  dt 
^2  c<t>i  dt 


+ 


502  502  dt  +  5^2^02  ^02  dt 

— C02  502  dt  4-  5^2502  C02  dt 

C02  C02  dt 


Let  [P]i  =  and  denote  the  averaged  version  of  these  forces  Since  the 

focus  here  is  on  the  longitudinal  equations,  only  and  are  of  interest.  Assuming  02  =  02  > 
these  average  forces  are 


^  =  ~  sin(^0/2)  sin 62  cos ^’2(^2  “  ^2) 

J=f  =  -^  sin(^0/2)  cos 6*2(02  -  02)- 


(3.16) 

(3.17) 


To  determine  the  aerodynamic  moment  due  to  the  force,  assume  that  [il2/i]i  =  [^2/1  2/2/1  so 
that 


^3/2- 


'  ^2/1  ' 

— C02502  -h  5025^2^02 

[-R3/i]i  — 

y2/i 

+ 

C02C02  +  5025^2502 

0 

502C^2 

On  the  down-stroke  the  aerodynamic  moment  is 


v/2 


msli  =  — ■7-pS3iy3/2<h)  C 


and  on  the  up)-stroke  the  moment  is 

V2 


msli  =  — T- P‘53 (2/3/202)  C 


C022/3/2  +  C022/2/1 
S022/3/2  -  2:2/1  C02 

-s02(s0?2/2/i +C0^a:2/i)  . 


C02  2/3/2  +  C022/2/1 
S022/3/2  -  a:2/lC02 

-502(502  2/2/1  +C02  2:2/1)  . 


(3.18) 


(3.19) 


Let  [m^Ji  =  (A^x,  Adx),  and  denote  the  averaged  version  of  these  forces  (Ad^,  Ady,  Ad®)*  For 
the  longitudinal  dynamics,  only  Aiy  is  required,  which  is 


Ml 


1 

T 


(2:2/1  sin(>l0/2)  -  2/3/2^0sin02)  cosP2(02  “  02)- 


(3.20) 


The  average  forces  and  moments  thus  far  calculated  are  for  a  single  wing.  When  the  wing  motion 
is  symmetric,  these  forces  and  moments  are  simply  doubled  in  magnitude.  That  is,  for  symmetric 
wing  motion,  the  total  average  forces  and  moments  are 


=  -^sin(2l^/2)sinP2COS02(02-02) 


2q 


sin(2l0/2)  cos  6*2(02  -  02) 

=  T  (^2/1  sin(2l0/2)  -  2/3/2^</>  sin  02)  cos  6*2(02  -  02)- 


(3.21) 

(3.22) 

(3.23) 
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or 


Tl  =  -2g/^ sin(>l^/2)  sin 02 cos ip2{^2  “  ^2) 

=  -2qf4,s\n{A^i2)  cos92{4>2  “  4>2) 

Ml  =  2g/^  {x2/i  sm{A^/2)  -  y3/2i4^sinV'2)  00862(^2  "  <;^2)- 
The  available  control  inputs  are:  02>  V’2>  02>  ^2- 


(3.24) 

(3.25) 

(3.26) 


3.4  Longitudinal  dynamics. 

Assuming  the  main  body  is  rigid,  then  with  the  standard  longitudinal  assumptions  {V  =  P  =  R  = 
0)  (Stevens  and  Lewis,  2003)  the  longitudinal  dynamics  are 

mil  =  —WQ  —  mgsinOi  +  Tx  (3.27) 

mW  =  UQ  +  mg  cos  0i  +  Tz  (3.28) 

lyQ  =  My.  (3.29) 

The  average  dynamics  axe  achieved  by  replacing  {Tx^  J^z'^My)  with  their  cyclic  averages  (.F°,  A1°). 

The  steady-state  dynamics  are  determined  by 

0  =  -W Q  -  mg sin  0i  -  2qf^ sm{A^/2)  sin  62  cos ip2{^2  “  ^2)  (3.30) 

0  =  UQ  +  mg  cos  0i  -  2qf^  sin(A,^/2)  cos  02(02  “  0“)  (3.31) 

0  =  “^.qf^,  (x2/i  sin(A,^/2)  -  y3/2A^sinip2)  cos 02(02  “  02)-  (3-32) 

For  any  given  steady  values  0i  ,U,W,Q  the  equilibrium  control  inputs  are 


.  _i /" 20:2/1  sin(A,^/2)>^ 


tan 


02 

$u;  =  - 


-1 


WQ  +  mg  sin  0i 


^—UQ  —  mg  cos  6i  cos  'ijj2 
—U Q  —  mg  cos  0i 


(3.33) 

(3.34) 

(3.35) 


2qf4,  sin{A^/2)  cos  62 ' 

To  determine  a  unique  set  of  control  inputs,  it  is  assumed  that  the  flapping  frequency,  f^,  and 
amplitude,  A^,  is  kept  constant.  From  this  specification,  it  follows  that 

(3.36) 


A(t,f(t>^2  ~  0202  "b  A^f^<p2  —  0. 

Since  both  02  ^nd  02  ^^e  positive  valued,  it  must  be  that  02  >  A^f^\  for  simplicity,  denote 

■■=  02  -  0^  (3.37) 

An  important  subset  of  steady-state  flight  is  the  hovering  mode  where  U  =  W  =  Q  =  0.  For  a 
given  pitch  angle,  ^i,  the  steady  actuation  requirements  are  given  by 

■  tan  0i ' 


Id,  lu 


62  =  tan 
02 


-1 


sm 


cos  t/^2 
2x2/1  sin(A,/,/2)  \ 

^01/3/2  / 

mg  cos^i 


2Q/,^sin(A,/,/2)cos02’ 


(3.38) 

(3.39) 

(3.40) 


where  apparently  02  and  02  must  be  non-zero. 
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3.5  Optimal  landing  trajectories. 


The  purpose  of  this  section  is  to  determine  a  trajectory  that  takes  the  vehicle  body  from  one  point 
to  another,  with  the  stipulation  that  the  ending  velocity  is  zero.  The  ending  point  is  considered 
fixed;  but  the  initial  conditions,  the  trajectory,  and  the  duration  are  all  open  for  analysis.  There 
are  obviously  infinite  number  of  trajectories  that  take  the  vehicle  from  point  A  to  point  B  while 
meeting  the  stipulated  requirements.  However,  generally  speaking,  this  should  be  accomplished 
with  as  little  ‘‘effort”  as  possible.  This  consideration  leads  to  optimization  problems,  the  solution 
of  which  will  specify  both  trajectories  and,  as  a  result,  the  forcing  requirements. 

Optimal  trajectories  are  derived  based  on  the  planar  inertial  equations 

Fx  =  mx  (3-41) 

Fz  —  mg  =  mz  (3.42) 

My  =  lyOl  (3.43) 

where  {x,-y,-z)  =  [ili]o,  (Fi,Fy,F^)  =  [f“]o,  and  {Mx,-My,-M^)  =  [m“]o.  The  forces  F^,  F^, 
and  moment,  My  are  inertial  representations.  The  conversions  are  [f“]i  =  Cio[f“]o  and  [m“]i  = 
Cio[m^]o.  Given  the  longitudinal  assumptions,  the  relationships  are  given  by 


^3.  =  Fx  cos  0  ~  Fz  sin  6 
Tz  =  ^Fx  sin  6  —  Fz  cos  6 
=  My. 


(3.44) 

(3.45) 

(3.46) 


The  objective  is  to  determine  the  trajectory  x(t),  y(t),  and  9{t),  subject  to  x{tf)  =  y{tf)  =  d{tf)  = 
0,  that  minimizes  a  specified  performance  index.  Extremum  of  a  prescribed  functional 

J  =  f  g(x,x,t)  dt 
Jo 

are  determined  by  a  solution  to  Euler’s  equation 

dt  \&k  J  5x 


Given  the  optimal  trajectory,  the  forces  and  moments  are  related  to  the  wing  kinematics  by 

'2x2/1  s\n{A4,/2)  2My  sin(>l,^/2)\ 


xj)2  =  sin 


(3.47) 


2qj^  sin(>l,^y 


(3.48) 

(3.49) 


where  ^  0. 

Since  the  forces  are  considered  independent,  the  minimization  of  Fx^  Fz^  My  can  be  performed 
independently.  The  functions  to  be  minimized  are 


(3.50) 

(3.51) 


(3.52) 


The  application  of  Euler’s  equation  results  in  2,  and  6  being  constant.  With  the  conditions 
x(0)  =  y(0)  =  ^(0)  =  0,  and  x{tf)  =  i(ty')  =  ^  =  0,  it  follows  that  the  optimal  trajectories  are 
given  by 


_  1 
Xf  4 

^(0  1 

Z;  4 

Of  4 


(3.53) 

(3.54) 

(3.55) 


where  Xf  :=  x(tf),  xq  :=  a:(0),  etc.  Each  of  the  trajectories  is  coupled  by  the  ending  time  per 


2xf  2zf  20  f 

io  io  60  ' 

Hence,  the  trajectories  can  also  be  written 

x{t)  =  — x/(r^  —  2r) 
x(t)  =  —Zf{T'^  —  2t) 

e{t)  =  -6»/(t^  -2t), 


(3.56) 


(3.57) 

(3.58) 

(3.59) 


where  r  :=  t/tf.  The  planar  trajectory  can  also  be  expressed  in  the  form 


z{x)  =  — X, 
Xo 


(3.60) 


which  indicates  a  linear  path.  The  path  is  either  ascending  or  descending  depending  on  the  sign  of 
io-  The  corresponding  forces  required  for  this  trajectory  are 


1 


r2 

Xq 


Fx  =  —7:^—  =  — 2m-f 

2  Xf  t) 


Xf 
1  Zq 

--m—  +  mg  ■■ 

Z  ^ 


cs 

2m-^  -h  mg 


t 


f 


M  -  -1/  ^  -  -27  ^ 

2  Hf~  Hy 


(3.61) 

(3.62) 

(3.63) 


The  required  inertial  forces  and  moment  are  thus  constant.  The  fixed  body  frame  forces  are  found 
by  Eqs.  (3.44-3.46). 

The  forces  and  moment  required  to  produce  the  optimal  trajectory  is  substantially  dependent  on 
the  ratios  xq/x/  and  z^jzf.  Recall,  that  an  important  approximation  in  deriving  the  aerodynamic 
forces  and  moment  is  that  |x(t)|  and  |i(t)|  are  “small”  relative  to  the  wing  motion.  The  wing 
kinematics  required  to  produce  these  forces  can  be  computed  by  Eqs.  (3.46-48). 


4  Landing  control. 

The  flight  control  problem  for  flapping  wing  MAV’s  has  been  addressed  in  several  works  [10].  Since 
the  flapping  motion  of  the  wings  The  consensus  approach.  This  section  details  control  methods 
that  are  suited  for  the  attachment  landing  problem.  These  techniques  are  applied  in  the  numerical 
example  of  Sec.  5.  Since  the  longitudinal  dynamics  are  significantly  nonlinear,  due  to  multiplicative 
and  trigonometric  terms,  the  focus  here  is  on  methods  of  nonlinear  control. 
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4.1  Control  model. 

To  construct  a  state-space  model  of  the  dynamics,  define  the  state  vector 


x=[x  z  e  U  W  Qf, 

and  the  control  vector 

u  =  [  02  V’2  ]^- 

Noting  the  inertial  position  dynamics 

X  =  U  cos  6  —  W  sin  0 
z  =  —U  sin  0  —  W  cos  0 
0i  =  -g, 

the  average  state  dynamics  (cf.  Sec.  4.2)  are  of  the  form 

Mt)  =  /(x)  +  g(u), 

where 

U  cos  0\  —  W  sin  9\ 
—U  sin^i  —  W  cosOi 

/(x)  = 

'  —WQ/m  —  gsinOi 

UQ/m  -f  5  cos  01 
0 


and 


(4.1) 


(4.2) 

(4.3) 

(4.4) 


(4.5) 


(4.6) 


g{u)  =[  0  0  0  J’l/m  J’j/m  My/Iy  ]  .  (4.7) 

From  Eqs.  (3.24-26),  the  forces  and  moment  are  related  to  the  control  input  by 

Tx  =  sin{A^/2)  sin  02  cos  ip2^w 

Tx  =  -2qf^  sm{A^/2)  sin  02  cos  ip2^w 
My  =  2qf^  [x2i\  sm{A^l2)  -  0.5>l,^j/3/2  sini/)2]  cB2^w 

To  produce  a  linear  form  of  the  model,  let 

Ax  =  X  -  xo  (4.8) 

Au  =  u  -  uo  (4.9) 

where  xo  is  an  equilibrium  point  of  Eq.  (4.5),  and  uo  is  the  corresponding  equilibrium  control 
input.  The  equilibrium  points  of  interest  are  those  of  hovering  flight,  i.e., 

Xo  =  [  0  0  0?  0  0  0  (4.10) 

UO=[02  V’2  (4.11) 
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The  linearized  dynamics  are  then 


Ax(t)  =  AAx(t)  +  BAu(t) 


where 

0 
0 
0 
0 

0 

and 


B  =  2qf^  sin(A^/2) 


-Usd^  -  Wcd^ 

cd1 

-sOl 

0 

-UcO^  +  Wse^i 

-c9° 

0 

0 

0 

0 

-1 

-gce^ 

0 

-Q/m 

-Wjm 

Qlm 

0 

0 

U/m 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0  ■ 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

644 

^45 

bie 

0 

0 

0 

bei 

bee 

bee 

0 

0 

0 

bei 

bee 

bee  . 

A  = 


where 


(4.12) 


(4.13) 


(4.14) 


644  =  -ce2Cip2^%/m 

645  =  se2Sip2^%/m 

646  = 

654  = 

655  =  0 

656  = 

bei  =  (0.5A^y3/2.s?/>2/sin(A^/2)  -X2/i)/4 
^65  =  -0.5A^y3/2c82Cip2^l,{sm{A^/2)Iy) 
bee  =  (0.5A^?/3/2Si/’2/sin(A^/2)  -  X2/i)ce^/Iy. 


The  linear  model  if,  of  course,  primarily  useful  when  the  states  are  kept  to  within  a  sufficient  mag¬ 
nitude  of  the  origin.  All  of  the  standard  techniques  of  linear  multivariable  control  are  then  available 
for  consideration.  For  the  present  case,  it  is  necessary  that  x,  z,  and  61  are  varied  over  a  large 
distance.  Extensional  linear  techniques  such  as  gain  scheduling  (Rugh,  1990),  and  simultaneous 
optimal  control  (Al-Sunni  and  Lewis,  1993),  are  available  in  these  circumstances.  Alternatively, 
under  certain  conditions  on  the  structure  of  the  equations  of  motion,  several  nonlinear  techniques 
are  available.  The  following  two  sections  apply  two  of  these  methods,  feedback  linearization  and 
backstepping,  to  the  flapping  control  problem. 


4.2  Feedback  linearization. 

Following  the  standard  approach  [cite],  define  the  output  vector 

y  =  [  X  2  9i  Y'  (4-15) 
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and  the  control  vector 


T={T^  T,  Myf  (4.16) 

Taking  the  second  derivative  of  the  output  results  in 

X  =  —W cos0(l  +  Q/m)  —  2gs[n6i  cos^i  —  U sin0i(l  -f  Q/m)  (4.17) 

-f  cos  QxTxjvn  —  sin^i^^/m 

z  =  W  sin0i(l  +  Qjvfi)  +  ^(sin^^i  —  cos^^i)  —  [/  cos^i(l  +  Qjm)  (4.18) 

~  sin  Q\Tx  j^n  —  cos  ^1.7^2 /m 

di  =  -Mylly.  (4T9) 

This  can  be  expressed  in  the  form 

y  W  =  /(y>  y)  +  Go(y,  y)F, 


where 


/(y,y) 


—W  cos0(l  +  Q/m)  —  2^  sin  cos0i  —  U  sin0i(l  +  Q/m) 
W sin0i(l  +  Q/m)  +5(sin^0i  —  cos^0i)  —  U  cos0i(l  +  Q/m) 

0 


(4.20) 


Go(y,y) 


cos  01 /m  —sin  01 /m  0 

—  sin0i/m  —  cos0i/m  0 

0  0  -1/4 


(4.21) 


Define  the  tracking  error 


e  =  y  -  r, 

where  r  is  the  reference  input.  The  force  input  is  then  defined  as 


dr  +  t{t) 


(4.22) 


(4.23) 


F  =  Go  1  -  Kie(t)  -  K2e(t)  +  K3 

where  Ki,  K2,  and  K3  are  diagonal  matrices  with  positive  elements.  The  error  dynamics  become 

e(t)  = '“Kie(t)  —  K2e(t)  —  K3  f  e{t)  dr,  (4.24) 

Jo 

which  are  clearly  exponentially  stable.  The  control  inputs  are  found  by 

u  =  gc(F)  (4.25) 

where  by  Eqs.  (3.47-49) 

1  f  2t2/i  sin(A^/2)  _  sin(A^/2)  \ 


(  ^tt>y3/2  } 


tan' 


-1  {  n  \ 

I  J’®COSl/>2  / 
2qf^  sin(A^/2)  cos  $2 


(4.26) 
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4.3  Backstepping. 

The  dynamics  are  divided  into  two  parts; 


17  =  fin)  + 

k  =  fain,0  +  GaF, 


where  77  =  (rr,  z,  ^1)  and  ^  =  ({/,  V,  Q);  thus 


fir,)  =  [  0  0  0 

fain)=[i-WQ/m-gsmei)  {UQ/m  + gcosOi)  0 


G  = 

Ga  = 


COS  0i  —  sin  0 
—  sin^i  —  cos^i  0 
0  0 

1/m  0  0 

0  1  /m  0 

0  0  1/4  ^ 


Define  the  function 


0  =  G-^ 


(4.27) 

(4.28) 


(4.29) 

(4.30) 

(4.31) 

(4.32) 


(4.33) 


where  Ki  and  K2  are  diagonal  matrices  with  positive  terms;  it  is  noted  that  G  ^  =  G.  Also,  define 
the  Lyapunov  function 

1/  =  l(a:2  +  +  ef).  (4.34) 

It  can  be  shown  (see  Sec.  14.4  of  [cite])  that  the  input 


F  =  G-i 


90 

dr] 


(/  +  G0)^ 


-fa-  m  -  <t>) 


(4.35) 


exponentially  stabilizes  the  system  of  Eqs.  (4.27,28).  As  before,  the  control  inputs  u  are  found 
from  Eqs.  (4.25,26). 


5  Numerical  example. 

The  preceding  sections  provide  the  analytical  framework  for  evaluating  the  attachment  landing 
requirements  and  methods  of  nonlinear  control  for  a  flapping  wing  MAV.  The  purpose  of  this 
section  is  to  demonstrate  the  application  of  these  results  with  a  numerical  example,  which  is  based 
on  the  a  hypothetical  vehicle  design,  resembling  the  general  planform  properties  of  the  bunblebee 
(Dudley  and  Ellington,  1990);  a  dynamic  flight  analysis  based  on  the  parameters  of  this  study  was 
conducted  by  Sun  and  Xiong  (2005).  The  baseline  design  parameters  are  given  in  Table  1. 
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Table  1.  Flapping  vehicle  parameters. 


Mass 

m  =  175  X  10  ^  kg 

Air  density 

p  =  1.25kg-m5 

Wing  reference  area 

5  =  53  X  10-5  j^2 

Moment  of  inertia 

ly  =  0.213  X  10-5  kg-m^ 

Flapping  amplitude 

=  116° 

Flapping  frequency 

=  155  Hz 

a:-distance  to  wing- root 

^2/1  =  3.91  mm 

distance  to  cop 

2/3/2  =  7.26  mm 

While,  as  previously  discussed  in  Sec.  4,  the  control  algorithms  are  based  on  the  average  dynam¬ 
ics,  subsequent  numerical  simulations  are  conducted  using  the  non-averaged  dynamic  model.  The 
basic  structure  of  the  simulation  model  is  shown  in  Fig.  2. 

5.1  Optimal  landing  trajectories. 

The  landing  requirements  are  calculated  based  on  the  optimal  trajectory  analysis  of  Sec.  3.5.  For 
this  example,  the  landing  point  is  taken  to  be  at  a;  =  0.1  and  z  —  0.05  m,  and  the  desired  pitch  angle 
is  =  —70°  (note  the  negative  angle  denote  a  pitch  angle  directed  toward  the  sky,  not  the  ground). 
Also,  the  initial  angular  rotation  rate  is  set  to  0^  =  0.01  deg/s.  The  resulting  actuation  requirements 
for  varying  landing  times,  tf,  are  shown  in  Fig.  3.  As  expected  the  actuation  requirements  increase 
as  the  landing  time  becomes  smaller.  It  is  important  to  note  that  these  trajectories  represent  open- 
loop  performance  based  on  the  averaged  dynamics.  Thus,  while  these  trajectories  are  theoretically 
achievable,  it  may  be  difficult  to  achieve  the  same  performance  in  feedback  control  for  the  non- 
averaged  system. 

5.2  Hovering  flight. 

One  of  the  primary  requirements  of  attachment  landing  is  the  ability  to  hover  at  a  large  pitch  angle. 
For  example,  consider  a  constant  hover  angle  of  6i  =  —70°,  in  which  case  the  A  and  B  matrices  of 
Eq.  (4.13,14)  are 


A  = 


B  = 


0 

0 

0 

-3.3552 

9.2184 

0 


0.3420 

0.9397 

0 

0 

0 

0 


0.9397 

-0.3420 

0 

0 

0 

0 


0 

0 

-1.0 

0 

0 

0 


0  0 
0  0 
0  0 
-2.9922  5.5561  x  10"3 
10.3369  0 

0  -2.1237  X  10^ 


0 

0 

0 

4.18879  X  10-5 
-3.9861  X  10-5 
0 


The  eigenvalues  of  the  A  matrix  are  all  zero.  This  is  the  case  for  any  hovering  equilibrium  (although 
inclusion  of  a  pitch  damping  term  results  in  a  negative  real  valued  eigenvalue).  Hence,  the  open-loop 
system  is  unstable,  and  must  be  stabilized  by  feedback. 
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The  demonstrate  the  the  basic  hovering  performance  with  classical  control  techniques,  a  standard 
linear  quadratic  regulator  optimal  control  (LQR)  is  applied  based  on  the  performance  index 


J  =  /  (Ax^QAx  +  Au^RAu)  dr. 

Jo 

After  some  iteration,  the  following  matrices  were  chosen: 


(5.1) 


Q  =  diag[l,  1,0.2, 0,0,0],  R  =  O.5I3. 


(5.2) 


Here  diag[-]  denotes  a  diagonal  matrix,  and  I3  is  a  3  x  3  identity  matrix.  The  hovering  performance 
is  shown  in  Figs.  4a, b  for  an  initial  pitch  angle  of  By  =  0.  The  results  show  that  the  linear 
controller  achieves  a  bounded  response.  The  performance  can  likely  improved  with  some  form  of 
integral  control.  However,  the  LQR  controller  was  found  insufficient  for  stabilizing  the  vehicle  at 
even  small  pitch  angles.  This  is  attributable  to  relatively  large  values  of  [/,  IT,  and  Q,  which  are 
assumed  small  in  the  linear  approximation. 


5.3  Landing  control. 

The  landing  sequence  is  divided  into  two  parts:  (1)  an  orientation  phase  where  the  pitch  angle  is 
brought  to  a  large  angle,  while  attempting  to  regulate  the  linear  motion;  and  (2)  a  terminal  landing 
phase  where  the  pitch  angle  is  kept  small  and  the  vehicle  is  brought  forward  to  its  target  landing 
position.  The  vehicle  is  initially  set  at  a  zero  pitch  angle  and  at  at  distance  of  100  mm  from  its 
landing  target.  The  desired  reference  pitch  angle  for  the  initial  orientation  phase  is  set  at  70°.  For 
the  following  simulations,  X2/1  is  set  to  zero.  This  was  found  to  produce  the  best  performance, 
likely  due  to  the  uncoupling  between  the  moment,  M.y,  and  the  vertical  force  Tz^ 

For  the  feedback  linearization  method  of  Sec.  4.3,  the  gain  matrices  of  Eq.  (4.24)  are  specified 
as 


Ki  =diag[100,100,4  x  10““*] 

K2  =  diag[20,20,0.4] 

K3  =  diag[10, 10,0.1]. 

These  were  chosen  based  on  numerous  iterations.  Generally,  it  was  found  that  the  desired  pitching 
rate  must  be  kept  relatively  small  to  maintain  a  reasonable  response  of  x  and  z  (within  10  mm 
of  the  origin).  Specifically,  for  the  particular  vehicle  under  consideration,  it  was  found  that  the 
inital  phase  must  be  in  excess  of  300  seconds.  Alternatively,  x  and  2  motions  can  be  tracked 
relatively  quickly.  The  numerical  simulation  of  the  the  two-part  landing  sequence  is  shown  Fig. 
5a;  the  corresponding  control  input  is  shown  in  Fig.  5b.  Once  a  pitch  angle  of  70°  is  achieved 
at  approximately  360  seconds,  a  forward  x  command  of  100  mm  is  initiated,  which  is  achieved  in 
approximately  5  seconds.  The  terminal  phase  of  the  landing  is  shown  in  Fig.  5c. 

For  the  backstepping  method  of  Sec.  4.3,  the  gain  matrices  of  Eq.  (4.32)  are  specified  as 

Ki  =  diag[100, 100,4  X  10"“*] 

K2  =diag[10,10,l]. 

The  two-part  landing  sequence  is  shown  in  Figs.  6a,  with  corresponding  control  inputs  shown  in 
Fig.  6b. 
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6  Conclusions. 


It  was  shown  that  three-degree-of  freedom  actuation  scheme  is  sufficient  for  the  attachment  landing 
problem.  These  degrees  of  freedoms  are  (i)  variation  of  of  down-stroke  and  up-stroke  flapping 
velocity;  (ii)  rotation  of  the  wing  at  its  root  attachment  point;  and  (iii)  wing  sweep.  The  rotation 
of  the  wing  is  especially  crucial  for  hovering  at  a  large  pitch  angle.  Also,  a  large  wing  rotation 
is  required  when  the  desired  landing  trajectory  is  “fast”.  Control  techniques  were  studied  using 
a  full  non- averaged,  nonlinear  numerical  simulation.  Methods  of  linear  control  were  found  to  be 
unsuitable  for  the  attachment  landing  problem;  a  demonstration  of  this  was  shown  with  a  standard 
LQR  method  for  stabilizing  hovering  dynamics.  The  nonlinear  methods  feedback  linearization  and 
backstepping  were  shown  to  be  suitable  for  landing  control.  The  application  of  these  methods  led 
to  a  two  part  landing  sequence,  the  first  part  of  which  is  a  large  pitch  rotation.  This  maneuver 
was  found  to  be  significantly  unstable,  and  could  not  be  performed  quickly.  The  second  part  of  the 
landing  sequence  is  a  forward  motion  (at  a  high  pitch  angle),  until  wall  contact  is  achieved.  This 
maneuver  was  found  to  be  significantly  more  stable,  and  could  be  performed  with  a  faster  settling 
time. 
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8  Appendix  -  Notation 

A  physical  vector  a  observed  from  two  reference  frames,  frame  with  orthonormal  unit  vectors 
(iojo^l^o)?  and  with  orthonormal  unit  vectors  can  be  expressed  as 

a  =  a^io  +  "b 


or 


a  —  cLx^i  +  “I" 

where,  in  general,  ^  Oy  ^  a'y^  ^  o!^.  We  will  refer  to  the  axes  of  frame  Sii  by  it  components 
Xu  yu  and  zi. 

The  matrix  representation  of  the  vector  a,  expressed  in  the  frame  is  denoted 


ab 


a^x 

ay 


Similarly,  the  frame  representation  of  a  is  denoted 


[a]i  = 


In  general,  a  symbol  representing  a  matrix  will  be  written  without  the  brackets.  However,  when 
it  is  necessary  to  identify  the  matrix  as  a  vector  representation,  we  will  use  the  bracket  notation 
[  ‘  ]i  where  i  denotes  the  corresponding  reference  frame.  In  general,  we  say  that  [a]i  is  the 
representation  of  a. 

The  matrix  [a]o  is  converted  into  the  matrix  [a]i  by  the  operation 


[a]i  =  Cio[b]o 


where  Cio  is  a  transformation  matrix  which  coverts  vector  representations  in  into  vector  rep¬ 
resentations  in  ^i.  The  form  of  this  matrix  is  dependent  on  the  specific  manner  in  which  rotation 
is  parameterized  (i.e.,  Euler  angles,  Euler  parameters,  etc.).  When  dealing  with  right-handed  or¬ 
thonormal  basis  vectors,  this  transformation  matrix  has  the  useful  property 

cio"  -  cTo 

That  is,  its  inverse  is  equal  to  its  transpose  (i.e.,  it  is  unitary). 

The  time  derivative  of  [a]i  is 

A([a]i)  =  A(C'io)[a]o  4-  CioA([a]o) 

where  the  time  derivative  is  denoted  Dt{-)  =  d/dt{*);  similarly,  we  write  At(*)  =  d^/dt^{-).  Pois¬ 
son’s  equation  provides  the  relation 


A(C'io)  =  ~[^i/o]iC'io 
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where  u^i/o  is  the  angular  velocity  vector  that  quantifies  the  rotational  motion  of  with  respect 
to  The  symbol  (*)  denotes  the  skew  symmetric  representation  of  the  vector.  That  is,  if 
[‘^i/oli  :=  [P  Q  P?'  then 


i/oli  — 


0  -R  Q 
R  0  -P 
-Q  P  0 


Since  Cio  is  an  orthonormal  operator,  we  also  have 

A(C'oi)  =  CoiliVifoll 

Thus,  the  time  derivative  becomes  of  [a]i  becomes 

A([a]i)  =  -[t^i/o]iC'io[a]o  +  CioDt{[a]o) 


or 


A([a]o)  =  Coi[u;i/o]I[a]i  +CoiA([a]i) 
The  second  derivative  follows  in  the  same  manner. 


9  Figures. 

Fig.  1.  Flapping  wing  kinematics. 

Fig.  2.  Simulation  model  schematic. 

Fig.  3.  Actuation  requirements. 

Fig.  4a.  LQR  hovering  control. 

Fig.  4b.  Actuation  inputs  for  LQR  hovering  control. 

Fig.  5a.  Feedback  linearization  control. 

Fig.  5b.  Actuation  requirements  for  feedback  linearization  control. 
Fig.  5c.  Terminal  landing  phase  for  feedback  linearization  control. 
Fig.  6a.  Backstepping  control. 

Fig.  6b.  Actuation  requirements  for  backstepping  control. 
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Fig.  1.  Flapping  wing  kinematics. 


Fig.  2.  Simulation  model  schematic. 
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Fig.  3.  Actuation  requirements. 


Fig.  4a.  LQR  hovering  control. 
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Fig.  4b.  Actuation  inputs  for  LQR  hovering  control. 
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Fig.  5a.  Feedback  linearization  control. 
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Fig.  5b.  Actuation  requirements  for  feedback  linearization  control. 
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Fig.  5c.  Terminal  landing  phase  for  feedback  linearization  control. 
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Fig.  6a.  Backstepping  control. 
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